% Расчёт КНД
% Цикл по фи(азимутальному углу)
Angls_phi=30; 
delta_phi=(pi/2)/Angls_phi;
phi_=-pi/2:delta_phi:0;
for ff=1:Angls_phi+1
        % РАСЧЁТ ДН РУПОРА
    clc
    M=10;  % Половинное число шагов интегрирование по осям
    Angls=50; delta=2*pi/Angls;    % Шагов по угловой координате и шаг для DN
    teta=-pi:delta:pi; % Полярный угол 
     %Матрица, соответствующая диаграмме направленности в плоскости (Z0PHI)
    Frup=ones(size(teta));
    %Расчет матрицы множителя системы
   [Frup]=d_raskrv(teta,phi_(ff),Angls, axr, byr, lmbd, M);
   maxFrup=max(abs(Frup));  Frup=abs(Frup./maxFrup);

    % Расчёт а2, а4
   n0_5=round(atan(tan(gama0)/2)/delta);
            f0_5=Frup(n0_5+Angls/2+1);

   n_1=round(gama0/delta);
            f_1=Frup(n_1+Angls/2+1);

   a4=(12+4*f_1-16*f0_5)/3;
   a2=f_1-1-a4;   % конец расчёта ДН рупора

   % Расчет ДН параболической антенны
   Angls=80;
   delta_p=pi/18/40;
   teta_p=-pi/18:delta_p:pi/18-delta_p; % Полярный угол параболической антенны
   F_parab=0;
   U=0;sin_teta=0;
   for kk=1:Angls/2
      U(kk)=(k*R0*sin(teta_p(kk)));
      L1=2/(U(kk))*besselj(1,abs(U(kk)));
      L2=8/(U(kk))^2*besselj(2,abs(U(kk)));
      L3=48/(U(kk))^3*besselj(3,abs(U(kk)));
      F_parab(kk)=((1+a2+a4)*L1-(a2/2+a4)*L2+(a4/3)*L3);
      sin_teta(kk)=sin(teta_p(kk));
   end
 F_p(ff)=max(F_parab);
    % Нормировка матрицы 
   maxF_parab=max(abs(F_parab));
   F_parab=abs(F_parab./maxF_parab);  % конец расчёта ДН параболической антенны
    % Расчёт интеграла по teta методом прямоугольных треугольников
   int_teta(ff)=sum(F_parab(1:kk).^2.*sin_teta(1:kk))*delta_p;
end               
Int=sum(int_teta)*delta_phi;
% Расчёт КНД
disp('КНД параболической антенны:')
D_max=abs(pi/(2*Int))
pause  
